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ABSTRACT 



A new approach to the stochastic analysis of general com- 



partment models is presented. The analysis is based on the concept 
of diffusion approximations. The state of a compartment system is 
represented as the superposition of a deterministic process, 
characterized by a system of ordinary differential equations, and 
a random noise process characterized by stochastic differential 
equations. All transition rate parameters are permitted to be time 
dependent. Numerical solutions are presented for the two-compartment 
case, and compared with those of Cardenas and Matis (1975a). 
Extensions to non-linear compartment models are discussed. 
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1. Introduction 



Quantitative representation of the distribution over time 
of a drug or pollutant in a human or animal body is now often 
carried out in terms of compartment models . Such elements as the 
blood, gut, liver, lean tissue, etc., are characterized as homo- 
geneous compartments, within which the drug resides for a time, 
later to transit to another compartment, perhaps recycling but 
eventually vanishing because of evacuation or metabolism mechanisms. 

Classical papers in this area, e.g. that by Bischoff et al. 
(1971) concerned with use of methotrexate in cancer therapy, pro- 
posed deterministic descriptions of flows between compartments, 
and stocks within, using differential equation systems. However, 
variations in drug concentrations over replicated experiments 
suggest a probabilistic or stochastic representation. Stochastic 
compartment analysis has undergone rapid development for several 
years, and might be categorized as follows. First, many papers 
have emphasized the mathematical formulation and exploration of the 
stochastic behavior of a variety of compartment models; papers of 
this type are exemplified by Bright (1973), Matis (1972, 1974, 

1975a, 1975b) , Marcus (1975) , Matis and Hartley (1972) , Purdue 
(1974a, 1974b, 1975) , Rescigno (1973) , Rubinow and Winzer (1971) , 
Thakur, Rescigno, and Schafer (1972, 1973), and Thron (1972). 

Second, several authors have investigated problems of statistical 
inference for the transition rates in stochastic compartment models; 
see Burkinshaw and Marshall (1971), Cornfield, Steinfeld, and 
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Greenhouse (1970) , Kodell and Matis (1976) , Matis and Hartley (1972) , 
Rodda et al. (1975), and Shah (1976). Last, the actual application 
of compartment models to biomedical, ecological, and chemical systems, 
as well as to pharmacokinetics, has been carried out, for example, 
by Metzler (1971), Rodenbeck et al. (1975), Sheppard (1962), Siegel, 
Cooper, and Meisner (1968) and Uppuluri and Bernard (1967). A large 
bibliography is given in Jacquez (1972). 

The present paper falls into the first category: we suggest 

a new approach to stochastic compartment analysis based on mathe- 
matical diffusion processes ; see Feller (1971) for an introduction, 
and Arnold (1973) and Gikhman and Skorohod (1971) for a systematic 
treatment. In short, we write stochastic differential equations 
to describe compartment concentration changes, and show that in a 
natural limit interesting joint distributions are jointly Gaussian 
or Normal. In the literature a variety of different compartment 
models have been formulated for one up to n-compartment , reversible 
or irreversible, open or closed, systems. These have generally 
been characterized by discrete-valued state variables for compart- 
ment contents, the latter changing according to multivariate 

birth-and-death processes. Kolmogorov forward equations for the 
transition probabilities are then derived, and solved by transforms. 

Purdue (1974) notes the similarity of these models to those 
for infinite server queues. We develop and discuss this connection 
in the appendix of this paper. For such discrete models it is possible 
to write down expressions for the moments of compartment contents at 
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a fixed time; under certain conditions (Poisson inputs, transition 
rates depend linearly on compartment contents) the contents of the 
compartments are statistically independent and Poisson distributed. 
Cardenas and Matis (1975a, b) have extensively investigated such 
models. The approach taken here, that of replacing a discrete birth 
and death process with a continuous diffusion, is approximate. How- 
ever, the results agree with the exact solution in the case of linear 
transition rates up to second moments. Furthermore, the limiting 
process may be shown rigorously to be normal, allowing the use of 
statistical inference procedures validated for normal distributions 
and processes. Perhaps more important is the fact that the diffusion 
approximation applies to situations with nonlinear transition 
Michaelis-Menten-type rates, where "exact" discrete-state methods 
yield cumbersome results; see McNeil and Schach (1973) for various 
examples concerning populations. It is reassuring to find that the 
diffusion approximations often agree exceptionally well with the 
birth and death solution; see Gaver and Lehoczky (1975,1976) for 
numerical comparisons. 



2. A Model 

We consider a general n-compartment model with time-dependent 
transition rates. Let C ± ( t) , i = 1 , . . . , n represent the contents 
of compartments l,...,n respectively at time t. We assume that 
the state of the system C(t) = (^ ( t) , . . . ,C R ( t) ) is a Markov process 
with transition probabilities given by 

P(A i (t)=l, A j ( t ) =0 , j ¥ i|C(t) ) = NX oi (t)dt + o(dt) 

i — 1 , . . . ,n 



P(A i (t)=-l, AjftJ-O, j * i|C(t>> = X tQ (t) C ± (t)dt + o(dt) (i) 

i = 1 / • • • t n 

P(A i (t)=-l, A j ( t ) =+l , A k (t)=0, k ? i,j|C(t)) = j ( t) C i (t)dt+ o(dt) 

for 1 £ if j < n, i ? j, A- > 0 . 
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where (t) = C^(t + dt) - C^(t), and N is a parameter that can 
be thought of as a dose level, eventually allowed to be large. 

The first of the three transition probabilities represents 
an increase in compartment i of size 1 from the outside. The 
second represents a decrease in compartment i of size 1 to the 
outside. Finally, the third represents a transfer of size 1 from 
compartment i to compartment j . 

The above formulation is rather general, and allows many 
different compartment systems to be studied simultaneously. It 
allows for an open system if Ag^(t) > 0 for some i = l,...,n 
or a closed system if AQ^(t) = ® f° r ^ = l' m ' n * The flows can 
be either reversible if A^j(t) > ® implies Aj^(t) >0 or 
irreversible if A^j (t) > 0 implies Aj^(t) = 0. 

We wish to study the case where (t) is large, say 

proportional to some number N. A diffusion approximation arises 
when N becomes large and we expand £(t) into terms of order N 
and /N. We can insure that C^(t) is large in two different 

ways. For open systems we assume, as given in (1), that the 
transition probabilities into the system from the outside are them- 
selves directly proportional to N. For closed systems where there 
is no input from the outside, we must assume that each compartment 
has contents proportional to N at time 0, that is C^(0) = Nc^(0). 

The effect of assuming C i (t) to be proportional to N 

is that transitions of all types occur very frequently. As a result, 
in any short time interval there will be many transitions of each 
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type. Each type of transition represents a change of +1 or -1 
in the contents of a certain compartment. In a short time interval 
the total change will be a sum of independent Bernoulli random 
variables, and, as such, normally distributed by virtue of the central 
limit theorem. Actually, much more is true. Focusing on the com- 
partment process {C(t) , t _> 0} over time rather than at a single 
fixed time we see that the process will have normally distributed 
increments and hence be Gaussian. This observation is the key to 
understanding the diffusion approximation approach. The mathematical 
development of this approach is not given here; see Kurtz (1971) 
and Barbour (1972) , and the paper by Schach and McNeil (1973) . 



3 . Mathematical Development 

In this section we derive the diffusion approximation approach 
outlined above in a manner that seems intuitively appealing. The 
exposition will be in terms of (Ito-type) stochastic differential 
equations; although, as will appear subsequently, the stochastic 
differential equations describing stochastic variations in compart- 
ment contents are not ambiguous of interpretation. The derivation 
to be presented next is supplemented by further discussion, reserved 
for the Appendix, that relates the present theory to the approach 
by Barbour (1974) , and also to infinite server queueing-type models. 
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Let, then, dC^(t) = C^(t + dt) - C^(t) represent the change 
in the contents of compartment i in time (t, t + dt) . The 
change dCL (t) may be viewed as a sum of independent random 
variables : 

(a) inputs from outside the system, with mean and variance 
N^(t) dt; 

(b) inputs from compartment j ( j ^ i) , with mean and variance 
X . . (t) C . ( t) dt; 

ji J 

(c) outputs to compartment k (k ^ i) , with mean and variance 
A ik (t) CL(t) dt; 

(d) outputs from the system via compartment i, with mean and 

variance X. rt (t) C. dt. 

lO 1 

Represent the change in contents of compartment i in the manner 
of diffusion, i.e. in terms of drift and diffusion terms: 



n 



n 



dC i (t) = [NX Qi (t) 



j=l 



I X . . (t) C. (t) + l X . . C . (t) ] dt 
_ n il 1 •_! D 1 D 



j=l 



n 




j=0 



( 2 ) 



n 



+ 7 v'A • • (t) C . (t) dW . . ( t ) 

£ 3 1 D ji v 



j=l 



i = 1,2, n; C (0) = N£(0) , and {VT ^ (t), 0 <_ i , j <_ n , i ^ j}. 



/ • • • r 



is a family of standard Wiener processes , see Arnold (1974) . 
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Expression (2) is motivated as follows: the drift term, in 

brackets, represents the expected change in (t) over (short) 
time interval (t, t + dt) , given the states of all compartments 
and the corresponding infinitesimal flow rates or transition 
probabilities. The remaining terms represent the various random 
components of change, modeled as Gaussian random variables; the 
latter is plausible as N -*■ °° on central limit theorem grounds. 

For instance, arrivals from outside the system in (t, t + dt) 

are Poissonian, and hence for large N, nearly Gaussian, 

with mean and variance specified by (a) above, and this accounts for 

the first bracketed and unbracketed terms in (2) . 

It is shown in Kurtz (1971) that 

C.(t) 

— - — -* c^(t) as N °° in probability; (3) 

see also Barbour (1974). By analogy with the central limit theorem, 

consider the normalized process 

C. (t) - N c. (t) 

X. (t) = — ; (4) 

/N 

hence express the concentration as 

C i (t) = Nc i (t) + /N X i (t) . (5) 

In matrix form, C(t) = Nc(t) + /N X(t) . 
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The vector function Nc(t) is referred to as the deterministic 

A-/ 

approximation, and /n X(t) is a stochastic noise process superimposed 

upon the deterministic approximation. It remains to characterize 

c(t) and X ( t) . The stochastic process X(t) = (X, (t) , . . . , X (t) ) 

/v rJ x n 

defined by (4) satisfies a stochastic differential equation similar 
to (2), and it can be derived from (2) by an application of Ito's 
Lemma (Arnold, p. 90, or Gikhman and Skohorod, p. 24). The stochastic 
differential equation governing X(t) is given by 



dX ( t) 



A ( t ) X(t) dt + B ( t) dW(t) 

nj +*> rJ 



- i/N (c 1 ( t) - 



(t) 



- A ( t) c ( t) ) dt ( 6 ) 

(V /N/ 



where 



A ( t ) = a (t) , 

rsj 1 J 



1 < i, j < n 



n 



with a . . ( t) = A . . ( t) for i ^ j , a. . (t) = - \ A . . (t) 

1 J J j = 0 



~0 ^ - (A Q1 (t),... , A Qn (t)) 



JW(t) — (WQ^(t)/.««/WQ^(t)/ W n Q ( t) / 2(t) t • • • / W 1 ^ ( t ) f 

' W n,n-1 (t) )T 



B(t) = (b ik (t)) 



1 < i < n, l<k< n(n+l) 



with 
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1 



1 



/ • • • / 



n 



b ii< t > - /X 0i TE > 



b i(„ + i) (t) “ - /A iO (t) C i (t,/N 



b i {2n+ (i-1) (n-1) +j-l) (t) b j (2n+(i-l) (n-l)+j-l) (t) 



= - /A i j (t) C i (t)/N , j > i 



b i (2n+ (i-1) (n-l)+j) (t) b j (2n+(i-l) (n-l)+j) (t) 



= - /A ij (t) C i (t)/N , 



j < i 



and all other b., =0. 

lk 

We now let N + °° in (6). First, the terms in B(t) 
involving /TTTTt) ( t) /N converge to (t) c^(t). Second, 

equation (6) contains a term proportional to N. The coefficient 
of this term must be identically 0 for all t >_ 0 or (6) will 
not converge to a finite limit. Consequently £(t) must satisfy 

c ' (t) = A n (t) + A ( t ) c(t) (7) 

ryj (J r*J 

and, given that c^t) satisfies (7) , equation (6) becomes 

dX ( t) = A ( t) X(t) + B ( t) dW ( t) (8) 

fsj /V f\J ryJ r>J 

where C^(t)/N terms are replaced by c^(t) in B^t) . Note that the 
modelling ambiguity sometimes associated with stochastic differential 
equation representations, see Arnold (1973), Chap. 10, does not arise 
because the diffusion coefficient, B(*), of (8) does not involve X. 
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The deterministic approximation to ^(t) is given by Nc(t) 
where c(t) is characterized by (7). This system of differential 
equations is stable in the sense of Liapunov because all eigenvalues 
of ^A(t) have strictly negative real parts. The stochastic noise 
process superimposed on the deterministic approximation is given by 
^Xtt) where X(t) is characterized by (8). The stochastic 
differential equation given by (8) describes a non stationary multi - 
variate Ornstein-Uhlenbeck process . Much is known about such 
processes, see Arnold (1974); we will return to a discussion of the 
noise process in a later section. 



4 . The Deterministic Approximation 

The system of differential equations given in (7) along 
with an initial configuration £(0) can be easily solved. If 
A.(t) = 0 , the closed system case, then 



t 

c(t) = exp (/ A(s)ds)"c(0) 



(9) 



2 

where exp(M) =I+M+M/2! + *** for any square matrix M. 
For the open system, nonhomogeneous , case with ^ ® one 
first multiplies by exp(- J A(s)ds) to obtain 



_d_ 

dt 



t 

(exp (-/ 

0 



A (s) ds) -c (t) ) 



t 

exp (-/ A (s) ds) (t) 



( 10 ) 
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This equation can be integrated and the result multiplied by 
t 

exp ( / A(s)ds) to give 
0 ~ 

t t t 

c(t) = / exp ( f A(s) ds) * A„ (u) du + exp ( / A(s)ds)*c(0) (11) 

~0u~ ~° 0" ^ 

Equation (11) provides a complete solution to the problem of 
describing c,(t) . Simplifications are not possible in general; 
however, equation (11) can always be solved numerically by standard 
numerical integration techniques. We can consider a special case 
which allows equation (11) to be substantially simplified. Suppose 
we assume that all transition probabilities, A^j (t) , 1 £ i £ n, 
0£j£n, i^j, exhibit the same time dependence, or that 
A. . (t) = A. .f(t) where f (t) is some function of time. The 
matrix A(t) now has the form 

A ( t ) = Af(t) (12) 

where 




Assume that A can be diagonalized and has eigenvalues 

0.,0 „,..., 0 , left eigenvectors Jl , and right eigenvectors 

1 z n ~n 

r. ,...,r . This gives 



ii 



A ( t ) = RD ( t) L with LR = I (13) 

fO ro/O ru A JfO rsJ ' ' 



where R= , £ 2 , . . . , r n ) , L= . . . , 1 ^) ' , and D(t) is a 

diagonal matrix with ith entry 0^f(t). Equation (13) can be 
integrated to give 



t 

/ A ( s) ds = RF(u,t)L 
u 



with EJ(u,t) a diagonal matrix having elements 
Finally (14) can be exponentiated to give 



t 

exp / A(s) ds = RG ( u , t ) L (15) 

u 



where £(u,t) is a diagonal matrix with the ith diagonal element 
t 

exp(0. / f(s)ds). Equation (15) can be plugged back into (11) 
u 

to give 



t 

c(t) = R, / S( u ' t )'fe^o (u)du + (16) 



We illustrate the use of (16) by considering an example 
discussed by Cardenas and Matis (1975a) . In that paper the authors 
consider a case with time-dependent transitions for a closed two- 
compartment system. For closed systems ^q(u) = 0 and only the 
second term of (16) need be calculated. In this case it is 
assumed that 



(14) 

t 

(0 i / f (s) ds) . 
u 



12 



X 01 (t) 



X 02 (t) - 0 



and 



*l 2 (t> 

A io 

A 21 (t) 



A 20 (t) 



= . 3/ ( t+1) , 
= .5/ (t+1) , 
= .4/ (t+1) 

f 

= . 6/ (t+1 ) 



A(s) = 



-.8 



.3 



-1 



s+1 



and 



A = 



-.8 



-1, 



f(s) = 



s+1 



The matrix A has eigenvalues 9, = -.9 + /. 13 , 

~ t 1 e. 

6- = -.9 - /. 13 . Also exp(9. / f(s)ds) = (1+t) 1 . 

1 0 

to calculate 



It is easy 



c(t) = 

aJ 



a/4 b/4 



(1+t) 



(1+t) 




•12/(.12+a ) .4a/(.12+a ) 

. 12/ ( . 12+b 2 ) . 4b/ ( . 12+b 2 ) 



■c (0) 



= -.1 + /. 13 



(17) 



b = -.1 - /7T3 



Equation (17) may be expanded to yield 
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c x (t) = ( . 6386750491c 1 (0) + . 5547001962c 2 (0) ) (1+t) 1 

6 2 

+ (. 361324909c 1 (0) - . 5547001962c 2 (0) ) (1+t) 

6 1 

c 2 ( t ) = ( .4160251472^ (0) + . 3613249509c 2 ( 0) ) (1+t) 1 

+ (-.4160251472c 1 (0) + . 6386750491c 2 (0) ) (1+t) 

with 

Q 1 = - .5394448725 
0 2 = -1.2605551275 



( 18 ) 



Transition probabilities for this model were approximated 
by Cardenas and Matis using a truncated infinite series approach. 
These transition probabilities can be derived directly from (18) 
by establishing appropriate initial conditions. For example, by 
setting c^(0) = 1, c 2 (0) = 0, c^(t) and c 2 (t) become p^(t) 
and P]_ 2 (t)f respectively. By setting c^(0) = 0, c 2 (0) = 1, 
c^(t) and c 2 (t) become p 2 ^(t) and p 22 (t), respectively. The 
results derived from (18) are compared with the approximations 
of Cardenas and Matis in Table 1. These results illustrate the 
accuracy of the deterministic approximations in (18) , since the 
two answers agree to at least 8 significant figures. Furthermore, 
the numbers derived from (18) always exceed the Cardenas and 
Matis results. This is reasonable since the latter approximation 
is derived by truncation and is hence an underestimate of the true 
value . 
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The results in (18) illustrate another very important 
point. Compartment modelling has been criticized because it only 
allows levels which are mixtures of exponentials in t when, 
in fact, these models often do not fit data well (See Marcus 
(1975) p. 339 and Matis (1972) • Much more complicated functions 
are required to give adequate fits including gamma, Pareto, 
Pareto-exponential, and first passage density functions. One can 
see that all of these types of functions can be produced by 
appropriate choice of f (t) . Mixtures of exponentials arise when 
f (t) = 1; however, an enormous class of level functions can be 
produced by varying choices of f (t) . The introduction of time 
dependent transitions not only adds realism but produces models 
which can fit real data sets. 
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p ll (t) 



p 22 (t) 



P21 (t) 



P 12 (t) 



t 


Equation (18) 


Cardenas and Matis 
4 -term approximation 


1 


.5902421829 


.5902421826 


2 


.4435620915 


.4435620872 


3 


.3652902915 


.3652902729 


4 


.3155676392 


.3155675836 


5 


.2807030964 


.2807029846 


1 


.5151767471 


.5151767467 


2 


.3596617846 


.3596617805 


3 


.2823115397 


.2823115195 


4 


.2356325933 


.2356325396 


5 


.2041834275 


.2041833201 


1 


.1501308716 


.1501308715 


2 


.1678006137 


.1678006133 


3 


.1659575075 


.1659575068 


4 


.1598700918 


.1598700879 


5 


.1530393378 


.1530393292 


1 


.1125981537 


.1125981539 


2 


.1258504603 


.1258504600. 


3 


.1244681306 


.1244681301 


4 


.1199025689 


.1199025659 


5 


.1147795033 


.1147794969 



TABLE 1: A Comparison of the Deterministic Approximation 

(18) with the Cardenas and Matis Four Term 
Approximation . 
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5. The Noise Process 



The stochastic elements of the compartment process (C(t), t >_ 0) 



are embodied in the nonstationary multivariate Ornstein-Uhlenbeck 
process described by (8) . We summarize a few key results about 
such processes. These results are well known and available in 
Arnold, Section 8.2. First, we assume X(0) = 0 that is the initial 

r<J pJ 

condition is nonstochastic and embodied in c(0). The marginal dis- 



that is an n-dimensional multivariate normal distribution where £ (t) 
is the unique nonnegative definite solution of the matrix Riccati 
equation 



tribution of X(t) is 




(19) 



E ' (t) = A ( t ) E (t) + E (t) A T (t) + B ( t ) B T (t) 

<v r>J ru A J /\j rJ r\J 



( 20 ) 



Since C(t) = Nc(t) + /N X(t) we find 




( 21 ) 



Let B ( t ) B T (t) = D ( t ) and D(t) = (d..(t)) with 

pj ro t\j l j 




( 22 ) 



17 



The covariance matrix E(t) can be calculated by various 
numerical methods applied to (20) . We illustrate the calculation 
in the case n = 2 studied by Cardenas and Matis . 



A(t) 



E(t) 



/ -U 10 (t)+X 12 (t)) X 21 (t) 

\ X l2^ - ^ X 20 ^ + X 21 

( On(t) °i2 (t) 
ai 2 (t) a 22 (t) 



(t) ) 



Equation (20) becomes 

s(t) = G ( t) s ( t ) + H ( t ) (23) 

AJ r\J nj /v 

where 



^s(t) = (a 1 ^(t) / a 12 (t), a 22 (t))' 



G(t) 

rJ 



'-2( X 1Q (t) + X 12 (t) ) 



X 1 2 ( t ) 



2X 21 (t) 

-(A 10 (t) + A 20 (t) + X 12 (t) + X 21 (t) ) 

2 X 12 (t) 



X 21 (t) 

-2( X 20 (t) + X 21 (t) ) 



B (t) 




+ 



+ 



( X 10 ( t ) + X 12^ t ^ C l^ t ^ + X 21^ t ^ C 2 
( X 12 ( t) c q ( t) + X 21 (t)c 2 (t)) 
Xl 2 (t) c i(t) + (x 2 Q (t) + X 21 (t)c 2 




Equation (23) can be solved in exactly the same manner as (7) with 
initial conditions s(0) = 0, that is C(0) is deterministic. We 



find 
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For the case considered by Cardenas and Matis 



s (t > - lit { 


1 -1.6 


! 

1-* 

• • 

00 00 


•- 


1 - lit s 




^ 0 


. 6 


-2 / 





The matrix G has eigenvalues <J>^, $3 where 



4>1 = 2d 1 = -1.978889745 



*2 = G 1 + 9 2 = - 1 * 8 ' 

<j>, = 2e, = -2.521110255 



k n ( u+ i) 



e i-! 



+ k 12 (u+1) 



V 1 



6-1 9-1 

H(u) =\ k 21 (u+l) + k 22 (u+l) 

e -1 0-1 

k 31 (u+l) L + k 32 (u+l) 



where the k. . coefficients are 
ID 

factor in the integrand is given 



computing from 
by 



(18) . 



The other 




with 



1 



1 



1 



R = 


[ .6513878188 


-.25 


-1.151387819 1 




\ .4243060905 


-.75 


1.32569391 / 



L 



.407957184 

.4615384615 

.1305557202 



.7085463502 

-.3076923077 

-.4998540425 



3076923077 \ 



-.6153846153 

1 






s(t) can now be easily found. Each of the three variance-covariance 
terms is a linear combination of terms of the form 

(j) . . — <j) . 

(1+t) 1 ((l+t) 3 1 - l)/(<f> j -0 i ) . 



The Ornstein Uhlenbeck process described by (8) also has a 
known covariance function k(s,t) = Cov(X(s), X(t) ) . The function 
k(s,t) is given by the equation (T denotes transpose) 



min (s , t) m m m 

k (s , t) = £(s) / $ (u) B (u) B (u) ($ (u) ) du $ (t) (25) 

/O r\J f\J r\J Aj 



with 



u 

$(u) = exp (/ A(s)ds) . 



We do not use the covariance function in this paper; however, 
this function is important in the development of a statistical 
analysis of data from a compartment model. One may have data, say 
readings of drug concentrations in the blood stream, and wish to 
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estimate the parameters of the model, the A^-'s. Given the data 
have a multivariate normal distribution, the likelihood function 
can be written and estimators derived. The papers of Hartley 
and Matis (1971) and Kodell and Matis (1976) provide analyses 
in special cases. The results of this paper indicate that these 
kinds of results can be carried much further into far more general 
compartment models. 



6 . Steady State 

In the case / Xg(t) ^ 0, an open system, the system will 
approach steady state if A^j (t) -*• A^j as t -*■ + °°. Steady state 
is a condition whereby the probability distribution which describes 
the marginal distribution of the system becomes time homogeneous 
even though the actual contents of the compartments continue to 
vary according to (8) . The steady-state solution for a simple example 
is described in the Appendix. 

The process {C(t), t >_ 0} in steady state can also be 
approximated by Nc + /N X(t) where 



0 = A^ + Ac or c 

r\J U s\J 



-A 

ro 




and {X(t), t > 0} satisfies 
— 



(26) 
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( 27 ) 



dX (t) = AX ( t) + B dW(t) 

rj njnj rJ r* 

a stationary multivariate Ornstein-Uhlenbeck process. This means 
that 

C(t) 'V (-NA _1 A , NE) (28) 

CU ' n ssj rJ U rJ 



where E is the unique nonnegative definition solution of 



T T 

AE + EA = -BB . 

rO rJ r\//\J 



(29) 



As an illustration we consider a two compartment open system, 
either reversible or irreversible 



A 

rJ 



BB 



T 



c 



(A 10 +A 12 ) 



21 



( A 20 +A 21^ 



A 01 + ( A 10 + A 12^ C 1 +A 21 C 2 



' (A 12 C 1 + A 21 C 2 ) 



A 02 + A 12 C 1 + ^ A 20 + A 21 ^ C 2 



A 20 +A 21 



l 12 



01 



l 21 



A 10 +A 12 / \ A 02 



^ A 20 +A 21 ^ ^ A 10 +A 12 ^ A 12 A 21 



(30) 



22 



0 



-1 



11 



2 (^]_o +A 12^ 



2 A 



21 



12 



22 



l 12 

0 



(X 10 +X 20 +A 12 +X 21 ) 



-A 



21 



BB' (31) 

rJnJ 



-2 A 



12 



2 ( A 20 +X 21* 



Equations (30) and (31) thus complete the description of 
the multivariate normal density given in (28) . While the n = 2 
case can be done in closed form, the general case can be easily 
solved by numerical methods. 



7 . Extension to Nonlinear Models 

The technique of diffusion approximations outlined and 
applied to general linear compartment models in the previous 
section also provides a tractable approach to handling nonlinear 
models. The compartment modelling literature is vast and the 
area is still in rapid development. Nevertheless, papers involving 
a stochastic process approach assume a linear system, that is, 
one in which transition rates from i to j are proportional to 
the contents of the ith compartment. It is unlikely that real 
pharmacokinetic processes operate so simply. Rather, the 
transition rates are more likely to be complicated functions of 
the contents of several compartment levels, and of time. Such 
complicated models have not appeared in the literature, perhaps 
because they present serious mathematical intractabilities using 
the standard Kolmogorov forward equations analysis. 
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The method of diffusion approximations can effectively 
handle nonlinear transition rates. The deterministic approximation 
will be a system of nonlinear ordinary differential equations. 

This will only rarely be solvable in closed form. Nevertheless, 
the numerical solution of such systems is a standard topic in 
numerical analysis. The important observation is that the stochastic 
differential equation representing the noise process superimposed 
upon the deterministic process will be nonstationary but linear 
(often Ornstein-Uhlenbeck) . This means that the diffusion approxi- 
mation approach will yield analytic results for nonlinear systems, 
since many results are readily available for linear stochastic 
differential equations (see Arnold, Chapter 8) . It is hoped that 
this approach will encourage pharmacokinetic researchers to con- 
sider introducing nonlinear models should an analysis of data 
indicate the need. 

Two additional points should be made. First, the analysis 
presented can be easily carried out when the process is represented 
in terms of volume and concentration rather than in terms of 
total contents. Second, the method of diffusion approximations 
illustrates that the compartment system contents will have a 
marginal Gaussian distribution. As such, one is in a position 
to give a statistical analysis of such a system either to estimate 
transition rates or, indeed, to design an experiment to make 
proper inferences about such parameters. 
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APPENDIX 



The purpose of this appendix is to further justify the 
earlier analysis, in particular, equations (7) and (8) and their 
implications. We do this in certain tractable simple cases, namely, 
those of one and two components. 



1 . The Theory of Kurtz (1971) and of Barbour (1974) . 

In Section 2 of Barbour (1974) , hereafter called B, 
assumptions are given that imply convergence of the sequence of 
processes 



i; N (t) = /N 



'S N (t> 

, N 



s (t) 



Sn 



(t) - Nc ( t) 

rJ 



/n 



N = 1,2,3,..., 



to the diffusion y(t), as N -*■ °° in D(0,T), i.e. the space of all 

rJ 

right-continuous functions with left-hand limits. Note that the 
y(t) of B will turn out to be the same as our stochastic noise 

rJ 

process X(t). We shall verify these assumptions for illustrative 

nj 

cases. 

First consider the sequence of continuous parameter Markov 
processes x N (t) = C N (t)/N on [0,T], where C N (t) is the many- 
compartment model whose transitions are described by our (1) . 

Notice that the changes in x N (t) are of size 0, + 1/N , and that 
the infinitesimal transition rates are of the form N ^gi' an< ^ 

N^ij ( C i (t)/N) , as specified by (2.1,B). We now discuss the assump- 
tion and the Theorem K of B for these illustrative situations. 



( A-l ) 
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Single-Compartment Model. Let (t) = A„, , and 

U X U X 

Aio(t) = A-^g constant and positive, and n = 1, which represents 
a single-compartment system. 

It is easily seen that there is a unique function 5 (t) 
satisfying, for 0 £ t £ T, the natural deterministic equation, 

( 2 . 2 , B) : 



5' (t) 




A 10 Ut) . 



The solution is, putting 5(0) = c(0), 



5(t) 



A 10*" 

c ( 0 ) e ■ LU + 



01 -X 10 t . 

( 1-e ) = 

10 



■( 



c ( 0 ) 



A 01 \ A 10 t A 01 

A 10 j A 10 ' 



which is unique; hence B's Assumption A is justified. Suppose 
0 < c(0) < Agi/ A io for example; other cases may be treated in 
analogous fashion. Then, in the terminology of Barbour (1974) , 



S 




c (0) 




(1 




and 



S £ = 



c (0 ) - e , c ( 0) 



-AT A 01 n A 10 T . 
e + t (1-e ) + e 

A 10 






where 0 < e < c(0) suffices. Clearly the transition rates are 
multinomial within S £ , and B's Assumption B is also justified. Now 
Theorem K follows, and implies that y N (t) = (x^t) - 5(t)) 

converges to the diffusion y(t), y(0) = 0, with the characteristic 

rJ ^ 



(A-2 



(A-3 



(A-4) 



(A-5) 
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iey (t) 



function <M0,t) = E[e 



] satisfying the differential equation 




(A-6 ) 



recognizable as describing the Ornstein-Uhlenbeck process 
whose stochastic differential equation is 



dy = -A 10 ydt + /A 01 + A 10 ?(t) dW ' 



( A-7 ) 



this being the same s.d.e. as that obeyed by our noise process 
X for the present setup. Hence the procedure leading to our 
equations (7) and (8) yields the same results as those rigorously 
established by Kurtz, adapted by Barbour. 

Note, too, that a steady— state solution will exist: simply 
set 5' = 0 in (A-2) to discover the deterministic component: 

C(°°) = A q 1 //A 10* The stochastic noise has, from (A-6) or (A-7), 
variance equal to the mean, namely, ^qi //A 10* t ^ us the steady-state 
compartment content is, according to our theory, approximately 



N becomes large. This is in accord 
with the fact that the exact distribution is Poisson ^NA Q1 /A 10 ) , 
as is well-known, and shown again subsequently in this Appendix. 
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Two-Component Model . Here C(t) - (C^(t), C 2 (t)), the 
allowed transitions and rates are summarized below. 

Transitions Transition Rates 

C^t), c 2 (t) + 



a) 


c 1 (t) + 1, 


O 

r+ 


O 

£ 


b) 


c 1 (t) - 1, 


c 2 (t) 


o 

t— 1 

r< 

£ 


c 1 (t) ^ 


N / 


c) 


c 1 (t) - 1, 


c 2 (t) + 1 


NXl2 ( 


Ci(t)\ 


N / 


d) 


c 1 (t) , 


C 2 (t) - 1 


NX 20 ( 


c 2 (t)\ 


N / 



State scaling, as in (2.1,B) alters the + 1 changes to + 1/N. 
The function £(t) = (£ (t) , (t) ) satisfies the differential 

equations 

^1 = X 01 “ ^10 + X 12 ^ ^1 

(A-9) 



V 2 X ± 2^1 X 20^2 



For simplicity, we have not allowed entries into compartment 2 
from outside, and have also omitted back flow from compartment 2 
to compartment 1 . The two equations are readily solved to produce 
a unique solution for all t > 0: 



q(t) = 



C 1 ( 0 ) 



exp [- (A 10 



x i2 ) t] 



l 01 



10 



+ ^ ex P ( < X io + x i2^ t ) 



12 
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^2 ^ = c 20 ex P(“ A 20 t ^ + 



01 



A 10 + X 12 



•j~- [1 “ exp(-X t) ] 
20 



C l (0) " X TT 
1 10 +A 12 



12 



A 20 “ (A 10 + A 12^ 



x [exp(-U 1Q + ^ 12 )t) " exp(-X 2Q t)] (A-10 ) 

Thus Assumption A is satisfied, with c^(0), 02 ( 0 ) > 0* Assumption B 
is also satisfied with 

0 < e < mint^m, I 2 (T)1 (A-ll) 

where 

L (T) = min £. (T) >0, i = 1,2; 

0 < t < T 



It is clear from the differential equations that such values will 
always exist. Hence 



&<*> ■ » 1/2 - *<*>) 

converges weakly to a bivariate diffusion. Proceeding further, 
evaluate the characteristic function of the limiting noise, 

E [exp (i6 1 y 1 (t) + i0 2 y 2 (t))] = $ ( 0 1 , 0 2 , t) ; (A-12) 

The latter satisfies the following partial differential equation, 
from (2 . 7 , B) : 
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[ e i A 10 + ( 0 1 + e 2 )A 12^ $ 0 1 9 2 A 20 $ e 2 



1 r„2 



2 t9 l A 01 + 9 l A 10^1 (t) + (_0 l +0 2 ) X 12 ^l (t)+ 9 2^20^2 ^ 



(A-13) 



This is recognizable as describing a bivariate Ornstein-Uhlenbeck 
process. Now one may differentiate (A-13) at 0^ = 0 2 = 0, 

utilizing the fact that 

V.(t) = E [y? (t) ] = -4> Q 0 (0,0, t), i = 1,2, 

x i i 

and 

V l 2 (t) = E[y 1 (t) y 2 (t)] = ~ $ e 1 e 2 ,0 



to derive differential equations for the elements of the variance- 
covariance matrix of y. Not surprisingly, these differential 
equations agree with (23) (with A^j constant, and, in this case, 

A 2 i = 0)/ i.e., V l( t) = c xl (t), V 12 (t) = a 12^ t ^ / V 2 = a 22 ^ * 

Thus if the same initial conditions are adopted for £ , c and 

iU ^ 

for and X the deterministic and noise components of the 

approximations of this paper and of B are seen to be identical. 
The validation can be extended with no difficulty in principle 
to the general situation. 
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2 . Relation to Infinite-Server and Linear Markov Population Models 



The class of compartment models described by the transition 
scheme (1) are similar to the classical Poisson-arrival infinite- 
server models of queueing theory and identical to certain Markov 
population models of Bartlett (1949 ) , and Kingman (1969 ) • We 
now briefly explore the connection, illustrating by the single- 
compartment and two-compartment examples. In particular, it will 
be shown that, even in the time-dependent parameter case, limiting 
Gaussian marginal distributions occur for compartment contents, 
in agreement with the development of the paper. 

Single-Compartment Model . Let C^(t,u) denote the number 
of arrivals to the system that have occurred in the time interval 
( u , t ) , 0 £ u £ t, and that are present in compartment 1 at time t. 

Let I (t,u) be the probability that an arrival occurring at 
time u is present in compartment 1 at time t. Consider the 
generating function of C^(t,u), 

C. (t,u) 

E [ ] = g(z- L ,t;u) (A-14) 

Now by the assumption of Poissonian arrivals, see (1) , we may 
write 

g(z^,t;u) = [z^I^(t,u)N\Q^(u) du + 1- I-^(t,u)NAQ^(u) du] g ( z^ , t ; u + du) , 

(A-15) 
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where the bracketted expression is the generating function of the 
contribution to C^(t,u) from arrivals during (u, u + du) ; 
independence permits the multiplcation . Subtract g(z^,t;u + du) 
from both sides, divide by du, and let du -*■ 0 . The differential 
equation 

dg (z. , t;u) 

- = g ( z^ , t; u) [ ( z-j^— 1) I 1 (t,u) NA Q1 (u)du] (A-16) 

results, the solution of which is, when u is set equal to zero, 

C, (t) t 

g(z 1 ,t;0)= E[z 1 |C 1 (0)=0] = exp[(z 1 ~l)N / ^ ( t ,u) X Q1 (u) du] , 

(A-17) 

the generating function of the Poisson distribution with mean 
t 

^ ^ XQ^(u)du. It now follows immediately from the central 

limit theorem that the distribution of 



x 1 (t) 



t 

c 1 (t) - N / I 1 (t,u) A Q1 (u)du 



/N 



(A-18) 



converges weakly to the normal law with mean zero and variance 
t 

/ I 1 (t,u) X Q1 (u)du. Furthermore, it is easy to verify that 



Var [y ( t) ] 



/ 

0 



I 1 ( t,u) A 0] _du = x 



01 



10 



(1 - e 10 ) 
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when X Q1 and X 1Q are both constants, and y(t) is the solution 
of the stochastic differential equation (A-7) which results from 
following the approximation program that is the topic of the paper; 
we start with C^(0) = 0. Verification in the time-dependent 
parameter case is also possible, but is more difficult. Thus the 
correct marginal distribution is veritable directly for this, 
the simplest, case. With more labor the joint distribution at 
different times may be found using an extension of the technique. 



Two-Compartment Model . Consider the generating function 
of (C^t/U), C 2 (t,u)), denoted by g ( z ± , z 2 , t ;u) . Then an argument 
analogous to that joint presented shows that g satisfies a 
differential equation, which when solved yields 

t 

g(z 1 ,z 2 ,t;0) = exp[N{(z 1 -l) / I 1 (t,u) X Q1 (u)du 

t 

+ (z,-l) / I 9 (t,u) X nl (u)du}] ( A-19 ) 

1 0 

implying that the contents of the two compartments are Poisson 
distributed at time t, and that C^(t) and C 2 (t) are independent. 
The same central limit theorem argument now shows joint normality 
as N *+• 00 in agreement with the claims of the body of the paper. 
Calculations omitted here show that the diffusion approximation pro- 
duces first and second moments that agree with those of the Poisson 
input model above. In particular, the correlation term, 0 2 , of (20) 
equals zero, signifying the independence that we noted in (A-19) . 
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